   This file documents the GNU Scientific Library.

   Copyright (C) 1996, 1997, 1998, 1999 The GSL Project.

   Permission is granted to make and distribute verbatim copies of this
manual provided the copyright notice and this permission notice are
preserved on all copies.

   Permission is granted to copy and distribute modified versions of
this manual under the conditions for verbatim copying, provided that
the entire resulting derived work is distributed under the terms of a
permission notice identical to this one.

   Permission is granted to copy and distribute translations of this
manual into another language, under the above conditions for modified
versions, except that this permission notice may be stated in a
translation approved by the Foundation.

Interface
=========

   All of the integration routines use the same interface.  There is an
allocator to allocate memory to hold control variables and workspace, a
routine to initialize those control variables, the integrator itself,
and of course a function to free the space when done.  For an
integrator algorithm, call it COOL (substitute any of the currently
existing algorithms) we then have

 - Function: gsl_monte_COOL_state* gsl_monte_COOL_alloc (size_t DIM)
 - Function: int gsl_monte_COOL_init (gsl_monte_COOL_state* S)
 - Function: int gsl_monte_COOL_free (gsl_monte_COOL_state* S),
 - Function: int gsl_monte_COOL_validate (gsl_monte_COOL_state* S,
     double *X_LOWER, double *X_UPPER,         unsigned long DIMENSION,
     unsigned long FUNCTION_CALLS)

 - Function: int gsl_monte_COOL_integrate (gsl_monte_COOL_state* S,
     double* X_LOWER, double* X_UPPER,         unsigned long DIMENSION,
     unsigned long FUNCTION_CALLS,         double* RESULT, double*
     ERROR,         ...)

   Notice the ellipses in the last argument to the actual integration
routine.  This is because the vegas algorithm (and perhaps others in
the future) returns extra information - in the case of vegas, it is the
\chi^2 of the result.

   In addition to the common function interface, the routines also share
the state variables

 - Control: gsl_rng* ranf
     which determines which random number generator will be used.
   and

 - Control: int verbose
     which says whether to print information about the calculation
     (though the actual use depends on the algorithm).
   In the near future, we expect that there will be a common interface
for selecting a "log" stream for reporting various information about
the performance of the particular algorithm.

   We describe the algorithm-specific variables.

   PLAIN: None

   MISER:

 - Control: double alpha
     alpha controls how the variances for the two sub-regions are
     combined.  The Numerical Recipes gang argue that for recursive
     sampling there is no reason to expect that the variance should
     scale like 1/N and so they allow the scaling to depend on \alpha

          Var(f) = {\sigma_a \over N_a^\alpha} + {\sigma_b \over N_b^\alpha}

 - Control: double dither
     Rather than exactly bisection the integration region, dither
     allows the user to introduce a bit of fuzz.  This helps in the
     case when the function to be integrated has some symmetry, say if
     it is peaked in the center of the hypercube.  If needed, dither
     should be around 0.1.

   VEGAS:

 - Control: double alpha
     For Vegas, `alpha' controls the stiffness of the rebinning
     algorithm: `alpha = 0' means never rebin.  It is typically set
     between one and two.

 - Control: double acc
     Setting `acc' allows vegas to terminate when the desired accuracy
     has been achieved, rather than after a certain number of function
     calls.  Setting it to a negative value disables this feature.

 - Control: long int max_it_num
     The maximum number of iterations to perform.

 - Control: int stage
     Setting this determines the "stage" of the calculation.  Normally,
     `stage = 0'.  Calling vegas with `stage = 1' retains the grid (but
     not the answers) from the previous run, so that one can "tune" the
     grid using a relatively small number of points and then do a large
     run with `stage = 1' on the optimized grid.  Setting `stage = 2'
     keeps the grid and the answers from the previous run and `stage =
     3' enters at the main loop, so that nothing is changed - this is
     like deciding to change `max_it_num' during the run.

 - Control: int mode
     The possible choices are `GSL_VEGAS_MODE_IMPORTANCE',
     `GSL_VEGAS_MODE_STRATIFIED',
     `GSL_VEGAS_MODE_IMPORTANCE_ONLY'.  This determines whether vegas
     will use importance sampling or stratified sampling, or whether it
     can pick on its own.  In low dimensions Vegas uses strict
     stratified sampling (more precisely, stratified sampling is chosen
     if there are fewer than 2 bins per box).

Example
=======

   Here we provide s simple example of using the integration routines.
Vegas was chosen at random.

     #include <math.h>
     #include <stdio.h>

     #include <gsl/gsl_math.h>
     #include <gsl/gsl_monte_vegas.h>

     double f0(double x[])
     {
       int i;
       double product = 1.0;

       for (i = 0; i < 10; i++)
         product *= x[i];

       return product;
     }

     main ()
     {
       double xl[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
       double xu[10] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1};

       double res = 0;
       double err = 0;
       double chisq = 0;
       int status = 0;
       unsigned long calls = 10000;
       unsigned long dimension = 10;

       gsl_monte_vegas_state* s = gsl_monte_vegas_alloc(10);
       gsl_monte_vegas_init(s);

       s->alpha = 1.5; /* default */
       s->verbose = 0; /* default, we prefer to remain ignorant! */
       s->stage = 0; /* default.  Start at the beginning */
       s->acc = -1.0; /* default. don't terminate when accuracy is reached */
       s->max_it_num = 5; /* default.  Do five iterations before quitting */

       status = gsl_monte_vegas_integrate(s, f0, xl, xu, dimension, calls,
                                         &res, &err, &chisq);
       if (status)
         printf("oops!\n");
       else
         printf("vegas(f0) = %f +- %f with \"chisq\" %f\n", res, err, chisq);

     }

The Future
==========

   In the future, the author of the Monte Carlo routines intends to add
more algorithms, greater control over the current ones,  more error
handling, and a more consistent interface.  Probably something like the
rng interface will evolve, so that something like

       integrator = gsl_monte_alloc(gsl_monte_vegas, dimension)

will be possible.  The old interface will probably stay around for a
while (though, since this is version < 1, the author does not want to
be held to such as statement).

Ordinary Differential Equations (ODEs)
**************************************

   This chapter describes functions for solving ordinary differential
equation (ODE) initial value problems.

IEEE floating-point arithmetic
******************************

   This chapter describes functions for examining the representation of
floating point numbers and controlling the floating point environment of
your program.  The functions described in this chapter are declared in
the header file `gsl_ieee_utils.h'.

Representation of floating point numbers
========================================

   The IEEE Standard for Binary Floating-Point Arithmetic defines binary
formats for single and double precision numbers.  Each number is
composed of three parts: a "sign bit" (s), an "exponent" (E) and a
"fraction" (f).  The numerical value of the combination (s,E,f) is
given by the following formula,

     (-1)^s (1.fffff...) 2^E

The sign bit is either zero or one.  The exponent ranges from a minimum
value E_min to a maximum value E_max depending on the precision.  The
exponent is converted to an unsigned number e, known as the "biased
exponent", for storage by adding a "bias" parameter, e = E + bias.  The
sequence fffff... represents the digits of the binary fraction f.  The
binary digits are stored in "normalized form", by adjusting the
exponent to give a leading digit of 1.  Since the leading digit is
always 1 for normalized numbers it is assumed implicitly and does not
have to be stored.  Numbers smaller than 2^(E_min) are be stored in
"denormalised form" with a leading zero,

     (-1)^s (0.fffff...) 2^(E_min)

This allows gradual underflow down to 2^(E_min - p) for p bits of
precision.  A zero is encoded with the special exponent of 2^(E_min -
1) and infinities with the exponent of 2^(E_max + 1).

The format for single precision numbers uses 32 bits divided in the
following way,

     seeeeeeeefffffffffffffffffffffff

     s = sign bit, 1 bit
     e = exponent, 8 bits  (E_min = -126, E_max = 127, bias = 127)
     f = fraction, 23 bits

The format for double precision numbers uses 64 bits divided in the
following way,

     seeeeeeeeeeeffffffffffffffffffffffffffffffffffffffffffffffffffff

     s = sign bit, 1 bit
     e = exponent, 11 bits  (E_min = -1022, E_max = 1023, bias = 1023)
     f = fraction, 52 bits

It is often useful to be able to investigate the behavior of a
calculation at the bit-level and the library provides functions for
printing the IEEE representations in a human-readable form.

 - Function: void gsl_ieee_fprintf_float (FILE * STREAM, const float *
          X)
 - Function: void gsl_ieee_fprintf_double (FILE * STREAM, const double
          * X)
     These functions output a formatted version of the IEEE
     floating-point number pointed to by X to the stream STREAM. A
     pointer is used to pass the number indirectly, to avoid any
     undesired promotion from `float' to `double'.  The output takes
     one of the following forms,

    `NaN'
          the Not-a-Number symbol

    `Inf, -Inf'
          positive or negative infinity

    `1.fffff...*2^E, -1.fffff...*2^E'
          a normalized floating point number

    `0.fffff...*2^E, -0.fffff...*2^E'
          a denormalized floating point number

    `0, -0'
          positive or negative zero

     The output can be used directly in GNU Emacs Calc mode by
     preceding it with `2#' to indicate binary.

 - Function: void gsl_ieee_printf_float (const float * X)
 - Function: void gsl_ieee_printf_double (const double * X)
     These functions output a formatted version of the IEEE
     floating-point number pointed to by X to the stream `stdout'.

   The following program demonstrates the use of the functions by
printing the single and double precision representations of the fraction
1/3.  For comparison the representation of the value promoted from
single to double precision is also printed.

     #include <stdio.h>
     #include <gsl/gsl_ieee_utils.h>

     main ()
     {
       float f = 1.0/3.0;
       double d = 1.0/3.0;

       double fd = f; /* promote from float to double */

       printf("   float ="); gsl_ieee_printf_float(&f); printf("\n");
       printf("promoted ="); gsl_ieee_printf_double(&fd); printf("\n");
       printf("  double ="); gsl_ieee_printf_double(&d); printf("\n");
     }

The binary representation of 1/3 is 0.01010101... .  The output below
shows that the IEEE format normalizes this fraction to give a leading
digit of 1,

        float = 1.01010101010101010101011*2^-2
     promoted = 1.0101010101010101010101100000000000000000000000000000*2^-2
       double = 1.0101010101010101010101010101010101010101010101010101*2^-2

The output also shows that a single-precision number is promoted to
double-precision by adding zeros in the binary representation.

Setting up your IEEE environment
================================

   The IEEE standard defines several "modes" for controlling the
behavior of floating point operations.  These modes specify the
important properties of computer arithmetic: the direction used for
rounding (e.g.  whether numbers should be rounded up, down or to the
nearest number), the rounding precision and how the program should
handle arithmetic exceptions, such as division by zero.

   Unfortunately there is no universal API for controlling these
features - each system has its own way of accessing them.  For example,
the Linux kernel provides the function `__setfpucw'
("set-fpu-control-word") to set IEEE modes, while HP-UX and Solaris use
the functions `fpsetround' and `fpsetmask'.  To help you write portable
programs GSL allows you to specify modes in a platform-independent
using the environment variable `GSL_IEEE_MODE'.  The library then takes
care of all the necessary machine-specific initializations for you when
you call the function `gsl_ieee_env_setup'.

 - Function: void gsl_ieee_env_setup ()
     This function reads the environment variable `GSL_IEEE_MODE' and
     attempts to set up the corresponding specified IEEE modes.  The
     environment variable should be a list of keywords, separated by
     commas, like this,

          `GSL_IEEE_MODE' = "KEYWORD,KEYWORD,..."

     where KEYWORD is one of the following mode-names,

          `single-precision'

          `double-precision'

          `extended-precision'

          `round-to-nearest'

          `round-down'

          `round-up'

          `round-to-zero'

          `mask-all'

          `mask-invalid'

          `mask-denormalized'

          `mask-division-by-zero'

          `mask-overflow'

          `mask-underflow'

          `trap-inexact'

          `trap-common'

     If `GSL_IEEE_MODE' is empty or undefined then the function returns
     immediately and no attempt is made to change the system's IEEE
     mode.  When the modes from `GSL_IEEE_MODE' are turned on the
     function prints a short message showing the new settings to remind
     you that the results of the program will be affected.

     If the requested modes are not supported by the platform being
     used then the function calls the error handler and returns an
     error code of `GSL_EUNSUP'.

     The following combination of modes is convenient for many purposes,

          GSL_IEEE_MODE=double-precision,mask-underflow,mask-denormalized

     This choice ignores any errors relating to small numbers (either
     denormalized, or underflowing to zero) but traps overflows,
     division by zero and invalid operations.

   To demonstrate the effects of different rounding modes consider
following the program computes e, the base of natural logarithms, by
summing a rapidly-decreasing series,

     e = 1 + 1/2! + 1/3! + 1/4! + ...
       = 2.71828182846...

     #include <math.h>
     #include <stdio.h>
     #include <gsl/gsl_ieee_utils.h>

     int main (void)
     {
       double x = 1, oldsum = 0, sum = 0;
       int i = 0;

       gsl_ieee_env_setup (); /* read GSL_IEEE_MODE */

       do
         {
           i++;

           oldsum = sum;
           sum += x;
           x = x / i;

           printf("i=%2d sum=%.18f error=%g\n",i, sum, sum - M_E);

           if (i > 30)
              break ;
         }
       while (sum != oldsum);

     }

Here are the results of running the program in `round-to-nearest' mode.
This is the IEEE default so it isn't really necessary to specify it
here,

     GSL_IEEE_MODE="round-to-nearest" ./a.out
     i= 1 sum=1.000000000000000000 error=-1.71828
     i= 2 sum=2.000000000000000000 error=-0.718282
     ....
     i=18 sum=2.718281828459045535 error=4.44089e-16
     i=19 sum=2.718281828459045535 error=4.44089e-16

After nineteen terms the sum converges to within 4 \times 10^{-16} of
the correct value.  If we now change the rounding mode to `round-down'
the final result is less accurate,

     GSL_IEEE_MODE="round-down" ./a.out
     i= 1 sum=1.000000000000000000 error=-1.71828
     ....
     i=19 sum=2.718281828459041094 error=-3.9968e-15

The result is about 4 \times 10^-15 below the correct value, an order
of magnitude worse than the result obtained in the `round-to-nearest'
mode.

   If we change to rounding mode to `round-up' then the series no
longer converges (the reason is that when we add each term to the sum
the final result is always rounded up.  This is guaranteed to increase
the sum by at least one tick on each iteration).  To avoid this problem
we would need to use a safer converge criterion, such as `while
(fabs(sum - oldsum) > epsilon)', with a suitably chosen value of
epsilon.

   Finally we can see the effect of computing the sum using
single-precision rounding, in the default `round-to-nearest' mode.  In
this case the program thinks it is still using double precision numbers
but the CPU rounds the result of each floating point operation to
single-precision accuracy.  This simulates the effect of writing the
program using single-precision `float' variables instead of `double'
variables.  The iteration stops after about half the number of
iterations and the final result is much less accurate,

     GSL_IEEE_MODE="single-precision" ./a.out
     ....
     i=12 sum=2.718281984329223633 error=1.5587e-07

with an error of O(10^-7), which corresponds to single precision
accuracy (about 1 part in 10^7).  Continuing the iterations further
does not decrease the error because all the subsequent results are
rounded to the same value.

References and Further Reading
==============================

   The reference for the IEEE standard is,

     ANSI/IEEE Std 754-1985, IEEE Standard for Binary Floating-Point
     Arithmetic

   A more pedagogical introduction to the standard can be found in the
paper "What Every Computer Scientist Should Know About Floating-Point
Arithmetic".

     David Goldberg: What Every Computer Scientist Should Know About
     Floating-Point Arithmetic. `ACM Computing Surveys', Vol. 23, No. 1
     (March 1991), pages 5-48

     Corrigendum: `ACM Computing Surveys', Vol. 23, No. 3 (September
     1991), page 413.

     See also the sections by B. A. Wichmann and Charles B. Dunham in
     Surveyor's Forum: "What Every Computer Scientist Should Know About
     Floating-Point Arithmetic". `ACM Computing Surveys', Vol. 24, No.
     3 (September 1992), page 319

Debugging Numerical Programs
****************************

This chapter describes some tips and tricks for debugging numerical
programs which use GSL.

Using gdb
=========

   Any errors reported by the library are routed through the function
`gsl_error'.  By running your programs under gdb and setting a
breakpoint in this function you can automatically catch any library
errors.  You can add a breakpoint for every session by putting

     break gsl_error

into your `.gdbinit' file in the directory where your program is
started.  If the breakpoint catches an error then you can use a
backtrace (`bt') to see the call-tree, and the arguments which possibly
caused the error.  By moving into the caller (`up', `up') you can
investigate the values of variable at that point.

   Here is an example from the program `fft/test_trap', which contains
the following line,

     status = gsl_fft_complex_wavetable_alloc (0, &complex_wavetable);

The function `gsl_fft_complex_wavetable_alloc' takes the length of an
FFT as its first argument.  When this line is executed an error will be
generated because the length of an FFT is not allowed to be zero.

   To debug this problem we start `gdb', using the file `.gdbinit' to
define a breakpoint in `gsl_error',

     bash$ gdb test_trap
     GDB is free software and you are welcome to distribute copies of it
      under certain conditions; type "show copying" to see the conditions.
     There is absolutely no warranty for GDB; type "show warranty" for details.
     GDB 4.16 (i586-debian-linux), Copyright 1996 Free Software Foundation, Inc...
     Breakpoint 1 at 0x8050b1e: file error.c, line 14.

When we run the program this breakpoint catches the error and shows the
reason for it.

     (gdb) run
     Starting program: /home/bjg/gsl/fft/test_trap

     Breakpoint 1, gsl_error (reason=0x8052b0d "length n must be positive integer",
         file=0x8052b04 "c_init.c", line=108, gsl_errno=1) at error.c:14
     14        if (gsl_error_handler)

The first argument of `gsl_error' is always a string describing the
error.  Now we can look at the backtrace to see what caused the problem,

     (gdb) bt
     #0  gsl_error (reason=0x8052b0d "length n must be positive integer",
         file=0x8052b04 "c_init.c", line=108, gsl_errno=1) at error.c:14
     #1  0x8049376 in gsl_fft_complex_wavetable_alloc (n=0, wavetable=0xbffff778)
         at c_init.c:108
     #2  0x8048a00 in main (argc=1, argv=0xbffff9bc) at test_trap.c:94
     #3  0x80488be in ___crt_dummy__ ()

We can see that the error was generated in the function
`gsl_fft_complex_wavetable_alloc' when it was called with an argument
of N=0.  The original call came from line 94 in the file `test_trap.c'.

   By moving up to the level of the original call we can find the line
that caused the error,

     (gdb) up
     #1  0x8049376 in gsl_fft_complex_wavetable_alloc (n=0, wavetable=0xbffff778)
         at c_init.c:108
     108           GSL_ERROR ("length n must be positive integer", GSL_EDOM);
     (gdb) up
     #2  0x8048a00 in main (argc=1, argv=0xbffff9bc) at test_trap.c:94
     94        status = gsl_fft_complex_wavetable_alloc (0, &complex_wavetable);

Thus we have found the line that caused the problem.  From this point we
could also print out the values of other variables such as
`complex_wavetable'.

GCC warning options for numerical programs
==========================================

   Writing reliable numerical programs in C requires great care.  The
following GCC warning options are recommended when compiling numerical
programs:

     gcc -ansi -pedantic -Werror -Wall -W -Wmissing-prototypes
       -Wstrict-prototypes -Wtraditional -Wconversion -Wshadow
       -Wpointer-arith -Wcast-qual -Wcast-align -Wwrite-strings
       -Waggregate-return -fshort-enums -fno-common -Wnested-externs
       -Dinline= -g -O4

For details of each option consult the manual `Using and porting GCC'.
The following table gives a brief explanation of what types of errors
these warnings catch.

`-ansi -pedantic'
     Use ANSI C, and reject any non-ANSI extensions.  These flags help
     in writing portable programs that will compile on other systems.

`-Werror'
     Consider warnings to be errors, so that compilation stops.  This
     prevents warnings from scrolling off the top of the screen and
     being lost.  You won't be able to compile the program until it is
     completely warning-free.

`-Wall'
     This turns on a set of warnings for common programming problems.
     You need `-Wall', but it is not enough on its own.

`-O4'
     Turn on optimization.  The warnings for uninitialized variables in
     `-Wall' rely on the optimizer to analyze the code.  If there is no
     optimization then the warnings aren't generated.

`-W'
     This turns on some extra warnings not included in `-Wall', such as
     missing return values and comparisons between signed and unsigned
     integers.

`-Wmissing-prototypes -Wstrict-prototypes'
     Warn if there are any missing or inconsistent prototypes.  Without
     prototypes it is harder to detect problems with incorrect
     arguments.

`-Wtraditional'
     This warns about certain constructs that behave differently in
     traditional and ANSI C. Whether the traditional or ANSI
     interpretation is used might be unpredictable on other compilers.

`-Wconversion'
     The main use of this option is to warn about conversions from
     signed to unsigned integers.  For example, `unsigned int x = -1'.
     If you need to perform such a conversion you can use an explicit
     cast.

`-Wshadow'
     This warns whenever a local variable shadows another local
     variable.  If two variables have the same name then it is a
     potential source of confusion.

`-Wpointer-arith -Wcast-qual -Wcast-align'
     These options warn if you try to do pointer arithmetic for types
     which don't have a size, such as `void', if you remove a `const'
     cast from a pointer, or if you cast a pointer to a type which has a
     different size, causing an invalid alignment.

`-Wwrite-strings'
     This option gives string constants a `const' qualifier so that it
     will be a compile-time error to attempt to overwrite them.

`-Waggregate-return'
     Warn if any functions that return structures or unions are defined
     or called.  Some older compilers might have problems with such a
     construct.

`-fshort-enums'
     This option makes the type of `enum' as short as possible.
     Normally this makes an `enum' different from an `int'.
     Consequently any attempts to assign a pointer-to-int to a
     pointer-to-enum will generate a cast-alignment warning.

`-fno-common'
     This option prevents global variables being simultaneously defined
     in different object files (you get an error at link time).  Such a
     variable should be defined in one file and referred to in other
     files with an `extern' declaration.

`-Wnested-externs'
     This warns if an `extern' declaration is encountered within an
     function.

`-Dinline='
     The `inline' keyword is not part of ANSI C. Thus if you want to use
     `-ansi' with a program which uses inline functions you can use this
     preprocessor definition to remove the `inline' keywords.

`-g'
     It always makes sense to put debugging symbols in the executable
     so that you can debug it using `gdb'.  The only effect of
     debugging symbols is to increase the size of the file, and you can
     use the `strip' command to remove them later if necessary.

Contributors to GSL
*******************

*Mark Galassi*
     Conceived GSL (with James Theiler) and wrote the design document.
     Wrote the simulated annealing package and the relevant chapter in
     the manual.

*James Theiler*
     Conceived GSL (with Mark Galassi).  Wrote the random number
     generators and the relevant chapter in this manual.

*Jim Davies*
     Wrote the statistical routines and the relevant chapter in this
     manual.

*Brian Gough*
     Wrote the FFT package and the relevant chapter in this manual.
     Also wrote the error handling infrastructure.

*Reid Priedhorsky*
     Wrote the root finding package and the relevant chapter in this
     manual.

*Gerard Jungman*
     Wrote the special function library, discrete hankel transforms,
     the rudimentary linear algebra functions, and the original version
     of the acceleration (Levin U) transforms.

*Mike Booth*
     Wrote the Monte Carlo library.

*Jorma Olavi Ta"htinen*
     Wrote the initial complex arithmetic functions.

*Thomas Walter*
     Wrote the initial heapsort routines and cholesky decomposition.

Copying
*******

   The subroutines and source code in the GNU Scientific Library
package are "free"; this means that everyone is free to use them and
free to redistribute them on a free basis.  The GNU Scientific
Library-related programs are not in the public domain; they are
copyrighted and there are restrictions on their distribution, but these
restrictions are designed to permit everything that a good cooperating
citizen would want to do.  What is not allowed is to try to prevent
others from further sharing any version of these programs that they
might get from you.

   Specifically, we want to make sure that you have the right to give
away copies of the programs that relate to GNU Scientific Library, that
you receive source code or else can get it if you want it, that you can
change these programs or use pieces of them in new free programs, and
that you know you can do these things.

   To make sure that everyone has such rights, we have to forbid you to
deprive anyone else of these rights.  For example, if you distribute
copies of the GNU Scientific Library-related code, you must give the
recipients all the rights that you have.  You must make sure that they,
too, receive or can get the source code.  And you must tell them their
rights.

   Also, for our own protection, we must make certain that everyone
finds out that there is no warranty for the programs that relate to GNU
Scientific Library.  If these programs are modified by someone else and
passed on, we want their recipients to know that what they have is not
what we distributed, so that any problems introduced by others will not
reflect on our reputation.

   The precise conditions of the licenses for the programs currently
being distributed that relate to GNU Scientific Library are found in
the General Public Licenses that accompany them.